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ABSTRACT 

We  have  applied  a  powerful  simulation 
methodology  known  as  ab  initio  crystal  prediction  to 
assess  the  ability  of  a  generalized  model  of  CHNO 
intermolecular  interactions  to  predict  accurately  crystal 
structures  and  densities  of  various  classes  of  explosives. 
174  crystals  whose  molecules  contain  functional  groups 
common  to  CE1NO  energetic  materials  were  subjected  to 
this  methodology  and  include  acyclic  and  cyclic 
nitramines,  nitrate  esters,  nitroaromatics,  and 
nitroaliphatic  systems.  The  results  of  these 
investigations  have  shown  that  for  148  of  the  174 
systems  studied  the  predicted  crystal  structures  matched 
the  experimental  configurations;  75%  of  these 
corresponded  to  the  global  energy  minimum  on  the 
potential  energy  surface.  The  success  rate  in  predicting 
crystals  with  structural  parameters  and  space  group 
symmetries  in  agreement  with  experiment  indicates  that 
this  method  and  interaction  potential  are  suitable  for  use 
in  crystal  predictions  of  similar  CHNO  systems  when  the 
molecular  configuration  is  known. 

1.  INTRODUCTION 

A  molecular  simulation  capability  of  particular 
importance  in  energetic  materials  development  is  that  of 
ab  initio  crystal  structure  prediction.  This  type  of 
calculation  has  the  potential  to  significantly  enhance  and 
optimize  the  design  and  development  of  advanced 
energetic  materials  because  it  allows  the  prediction  of  a 
key  property  associated  with  performance  of  the 
material:  its  crystal  density.  The  ability  to  predict  the 
crystal  density  will  allow  a  designer  to  assess  detonation 
or  gun  performance  without  the  need  to  synthesize  the 
material.  In  this  manner,  candidates  with  densities  that 
will  produce  less-than-satisfactory  performances  can  be 
eliminated  from  further  consideration  before  actual 
resources  are  expended  on  synthesis.  Therefore,  this 
important  methodology  is  considered  to  be  a  key 
computational  tool  in  the  development  of  advanced 
energetic  materials  research. 


In  addition  to  allowing  a  materials  designer  to  rapidly 
examine  a  wide  variety  of  notional  candidate  materials  for 
applications  in  which  crystal  morphology  or  density  is  a 
critical  indicator  of  behavior,  activity  and  performance,  the 
method  can  be  used  to  assess  molecular  simulation  models 
that  might  be  used  in  other  types  of  molecular  simulations 
(i.e.  molecular  dynamics). 

The  success  of  the  method  is  dependent  on  the  quality 
of  the  model  used  to  describe  intermolecular  interactions 
among  the  atoms  in  the  crystals.  The  model  used  in  this 
study  is  a  generalized  CHNO  intermolecular  interaction 
potential  developed  for  use  in  other  types  of  molecular 
simulations  (Sorescu  et  ah,  1997).  This  model,  denoted  as 
the  SRT  model  after  its  original  authors,  has  already  been 
subjected  to  a  series  of  molecular  simulation  studies  in 
which  its  ability  to  predict  properties  of  a  variety  of 
molecular  crystals  has  been  evaluated  (Sorescu  et  ah, 
1998a,  1998b,  1998c,  1999a,  1999b).  Although  the  SRT 
model  appears  to  be  a  promising  model  for  use  in 
molecular  simulation  of  CHNO  molecular  crystals, 
increasingly  rigorous  assessments  should  be  performed  to 
establish  the  limits  of  the  model’s  capabilities  and  provide 
information  for  refinement.  Particularly,  it  is  important  to 
establish  whether  such  a  potential  is  able  to  reproduce  the 
crystal  structure  and  the  corresponding  lattice  energy 
known  from  experimental  data  and  also  to  predict  the 
crystallographic  symmetiy  without  its  prior  knowledge. 
The  study  presented  here  provides  such  an  evaluation.  In 
this  work,  we  have  subjected  a  series  of  CHNO  crystals 
that  are  representative  of  different  classes  of  explosives  to 
ab  initio  crystal  structure  prediction  using  the  SRT 
interaction  potential,  in  order  to  appraise  the  performance 
of  this  method  and  model  in  generating  hypothetical 
crystal  structures  that  are  consistent  with  those  observed 
experimentally.  We  have  attempted  to  make  this  analysis 
comprehensive  by  considering  a  large  number  of  chemical 
systems  of  interest  for  energetic  materials  applications.  In 
particular,  in  this  work  we  have  considered  compounds 
containing  nitramines,  nitrate  esters,  nitroaliphatics, 
nitroaromatics,  and  furoxans. 
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2.  EXPERIMENTAL  DATA  USED  IN  THE 
STUDY 

The  structural  data  of  the  molecular  systems 
considered  in  this  study  were  obtained  from  the 
Cambridge  Structural  Database  (CSD)  (Allen,  2002).  No 
modifications  to  the  experimental  structures  were  made. 
The  main  categories  of  compounds  used  were  species 
having  nitro  groups  attached  to  carbon  (denoted  as 
nitroaliphatic  or  nitroaromatic),  nitrogen  (nitramine),  or 
oxygen  (nitrate  esters).  Additionally,  a  few  furoxan 
systems  were  considered.  Within  this  paper,  systems 
will  be  referred  to  by  their  CSD  entry  identifier 
{refcode).  Structures  were  eliminated  from  consideration 
for  this  study  if:  1)  The  space  group  could  not  be  treated 
by  the  software  used  in  the  calculations,  MOLPAK 
(Holden  et  al.,  1993;  Busing,  1981;  Cromer  et  al.,  1987). 
We  did  not  exclude  crystals  with  space  groups  that  were 
alternate  settings  of  those  that  could  be  treated  by  the 
version  of  MOLPAK  used  in  this  study.  2)  The  systems 
contained  atoms  other  than  C,  H,  N  and  O.  3)  The 
systems  contained  solvent  molecules.  4)  There  is  more 
than  one  molecule  in  the  asymmetric  unit.  5)  Those 
systems  in  which  the  asymmetric  unit  contains  only  part 
of  the  molecule  and  no  C2  rotation  axis  or  center  of 
inversion  exist  among  the  symmetry  elements  of  the 
crystal.  However,  when  such  symmetry  elements  are 
present  the  corresponding  crystals  can  be  considered  for 
MOLPAK  simulations.  5)  Fractionals  for  some  or  all  of 
the  non-H  atoms  were  not  available.  6)  Systems  that 
have  disordered  atomic  positions  or  have 
crystallographic  R- factors  greater  than  1 1%. 

This  search,  which  was  not  exhaustive,  resulted  in 
the  identification  of  174  molecular  crystals  that  represent 
our  working  database  to  be  subjected  to  ab  initio  crystal 
prediction.  Information  about  these,  including  full 
chemical  name,  CSD  refcode,  chemical  formula,  space 
group,  temperature  at  which  measurements  were  taken, 
and  reported  R-factors  are  given  in  Rice  and  Sorescu 
(Rice  and  Sorescu,  2004). 

3.  INTERMOLECULAR  POTENTIAL 

The  PES  is  represented  by  a  very  simple  model  for 
use  in  rigid-molecule  simulations;  its  form  is  a  simple 
atom-atom  exponential-six  plus  Coulombic  potential. 
The  partial  charges  used  in  the  Coulombic  interaction 
terms  were  determined  by  fitting  atom-centered  partial 
charges  to  a  quantum-mechanically  determined 
electrostatic  potential  for  a  single  molecule  whose 
structure  corresponded  to  that  in  the  crystal  at  ambient 
conditions.  The  remaining  exponential-six  parameters 
were  adjusted  to  reproduce  the  experimental  structure  of 
the  RDX  crystal  at  ambient  conditions  (Sorescu  et  al., 
1997). 


4.  COMPUTATIONAL  METHOD 

4.1.  The  MOLPAK/WMIN  Crystal  Structure 
Prediction  Method 

In  our  study,  we  have  used  the  ab  initio  crystal 
structure  prediction  software  MOLPAK/WMIN  (Holden 
et  al.,  1993;  Busing,  1981;  Cromer  et  al.,  1987),  developed 
with  an  emphasis  on  density  predictions  of  energetic 
materials.  A  key  assumption  in  the  MOLPAK/WMIN 
procedure  is  that  of  closest  packing  of  the  molecules  in  the 
crystal;  the  procedure  attempts  to  identify  structures  with 
the  highest  densities.  The  MOLPAK  portion  of  the 
calculation  generates  numerous  candidate  crystal  structures 
of  minimum  volume  for  molecular  coordination 
geometries  observed  in  the  most  common  space  groups. 
MOLPAK  considers  29  coordination  geometries  for  the 
triclinic  (PI,  P\  ),  monoclinic  (P2b  P2!/c,  Cc,  C2,  C2/c) 
and  orthorhombic  (Z=4,  P2i2]2,  P2j2i2i,  Pca2i,  Pna2[  and 
Z=8,  Pbcn,  Pbca)  space  groups.  A  subset  consisting  of  the 
most  dense  candidate  crystals  are  subjected  to  further 
refinement  through  full  energy  minimization  using  the 
program  WMIN  (Busing,  1981).  The  steps  in  the 
procedure  are: 

1.  The  Cartesian  coordinates  of  the  three- 
dimensional  molecular  model  that  is  to  be  used  in  building 
the  candidate  crystals  (hereafter  referred  to  as  the  “test 
molecule”)  are  specified.  The  centroid  of  the  test  molecule 
is  located  at  the  origin  of  the  hypothetical  crystal.  The 
orientation  of  the  molecule  is  given  by  a  set  of  Euler 
angles.  In  this  study,  the  arrangement  of  the  atoms  in  the 
molecule  corresponds  to  that  of  the  crystal  molecular 
structure,  as  reported  in  the  CSD  (Allen,  2002).  An 
alternate  method  of  determining  the  coordinates  of  the  test 
molecule  is  through  quantum  mechanical  prediction. 

2.  An  initial  packing  arrangement  is  obtained  by 
surrounding  the  test  molecule  with  a  coordination  sphere 
containing  other  molecules.  The  contents  and  three- 
dimensional  structure  of  the  sphere  are  dependent  on  the 
crystalline  space  group  symmetry.  The  definitions  of  the 
various  coordination  spheres  used  in  a  MOLPAK 
calculation  were  obtained  from  detailed  analyses  for  a 
large  number  of  organic  crystal  structures  (Holden  et  al., 
1993).  The  analyses  showed  that  the  most  probable 
number  of  molecules  in  the  coordination  sphere  is  14,  and 
that  specific  “patterns  and  sub-patterns”  were  apparent  in 
the  three-dimensional  structure  of  the  coordination  spheres 
(Holden  et  al.,  1993). 

3.  Adjacent  molecules  in  the  coordination  sphere  are 
then  systematically  moved  in  small  steps  toward  the 
centrally-located  test  molecule.  At  each  step,  the  potential 
energy  of  the  “crystal”  is  calculated.  This  continues  until  a 
repulsion  criterion  is  met.  Once  this  criterion  is  met,  the 
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packing  procedure  stops,  and  the  volume, 
crystallographic  parameters  and  Eulerian  angles  that 
describe  the  orientation  of  the  test  molecule  are  stored. 

4.  The  test  molecule  is  then  reoriented  about  only 
one  of  the  Eulerian  axes  by  10°,  and  the  packing 
procedure  (steps  1-3)  is  repeated  until  the  entire  Eulerian 
space  is  sampled.  Rotations  in  10°  steps  about  the  three 
Eulerian  axes  will  result  in  the  generation  of  6,8  5  9  (193) 
orientations  and  hypothetical  crystal  structures  for  each 
of  the  possible  space  group/coordination  sphere 
combinations. 

5.  After  the  -7000  structures  were  generated  for 
each  coordination  sphere  geometry,  they  are  ranked 
according  to  density.  The  25  most  dense  structures  are 
subjected  to  full  energy  minimization,  where 
minimization  is  performed  with  respect  to  the 
crystallographic  parameters.  The  energy  minimization  is 
space-group  symmetry  restricted  (i.e.,  the  space  group 
symmetry  is  conserved  throughout  the  minimization), 
and  performed  using  the  code  WMIN  (Busing,  1981). 

At  the  completion  of  the  energy  minimizations  of 
the  25  most  dense  hypothetical  crystals  generated  for 
each  of  the  29  space  group/coordination-sphere 
geometries,  the  user  has  the  information  needed  to 
identify  the  “correct”  crystal  structure  according  to  his 
specific  criteria  (e.g.  lowest  lattice  energy  or  highest 
density  [Lommerse  et  al.,  2000]).  During  energy 
minimizations  some  of  the  hypothetical  crystals  can 
converge  to  the  same  local  minimum  on  the  PES. 
Elowever,  identical  crystals  have  not  been  eliminated  in 
the  analyses  reported  hereafter. 

The  MOLPAK/WMIN  procedure  proceeds  under 
the  assumption  that  all  important  regions  of 
configuration  space  are  adequately  sampled  and  that  the 
interatomic  interactions  are  sufficiently  accurate  to 
describe  packing  and  lattice  energies.  Unfortunately,  we 
found  that  in  some  cases,  the  MOLPAK/WMIN 
procedure  did  not  result  in  the  identification  of  the 
experimental  structure.  For  those  cases,  a  second 
MOLPAK/WMIN  series  of  calculations  were  performed 
in  which  the  test  molecule  was  reoriented  about  each  of 
the  Eulerian  axes  in  increments  of  5°,  resulting  in  the 
generation  of  50,653  (37')  orientations  and  a 

corresponding  larger  number  of  hypothetical  crystal 
structures.  For  a  number  of  structures,  the  increased 
resolution  of  the  grid  led  to  successful  identification  of 
the  corresponding  experimental  structures. 

4.2.  METHOD  OF  COMPARISON  OF 
STRUCTURES 

The  comparison  of  crystallographic  structures 
predicted  using  the  MOLPAK/WMIN  codes  with 


experiment  was  accomplished  in  a  series  of  steps.  The 
first  step  is  the  comparison  of  the  reduced  cell  parameters 
obtained  from  MOLPAK  calculations  with  the 
corresponding  experimental  values.  The  MOLPAK  crystal 
with  reduced  cell  parameters  that  most  closely  agree  with 
experimental  values  is  identified  as  a  potential  match. 

A  more  rigorous  comparison  of  crystals  has  been 
performed  to  check  that  not  only  the  lattice  parameters  are 
in  agreement  with  experimental  values  but  also  that  the 
arrangement  and  orientation  of  molecular  systems  inside 
the  unit  cell  resemble  the  experimental  situation.  The  first 
step  of  this  portion  of  the  assessment  is  the  generation  of 
the  experimental  unit  cell  and  a  MOLPAK  supercell 
consisting  of  125  unit  cells  (5x5x5  unit  cells).  Next,  the 
experimental  unit  cell  and  MOLPAK  supercell  are 
translated  such  that  the  mass  centers  of  the  central 
molecule  of  the  MOLPAK  supercell  and  one  of  the 
molecules  in  the  experimental  unit  cell  are  located  at  the 
origin.  Using  the  procedure  described  by  Kearsley 
(Kearsley,  1989),  the  MOLPAK  supercell  is  then  rotated 
about  the  origin  such  that  the  central  molecule  of  the 
MOLPAK  supercell  is  superimposed  onto  the 
experimental  molecule.  Since  the  molecular  structures 
used  in  the  MOLPAK  calculations  are  exactly  the  same  as 
those  in  the  experimental  unit  cells,  the  Kearsley  procedure 
will  result  in  a  perfect  superposition  of  the  MOLPAK  and 
experimental  molecules.  At  this  point,  the  molecules  in 
the  MOLPAK  supercell  that  correspond  to  each  of  the 
remaining  equivalent  molecules  in  the  experimental  unit 
cell  are  identified  by  minimizing  the  root-mean-square 
(rms)  deviation  of  distances  between  the  two  sets  of 
molecules.  Once  the  MOLPAK  unit  cell  that  most  closely 
resembles  the  experimental  unit  cell  is  obtained,  the 
MOLPAK  and  experimental  unit  cells  are  shifted  such  that 
the  mass  centers  of  both  cells  are  located  at  the  origin.  The 
Kearsley  procedure  (Kearsley,  1989)  is  then  applied  using 
all  molecules  within  both  unit  cells,  to  obtain  optimum 
superposition  of  the  MOLPAK  unit  cell  onto  the 
experimental  cell.  Locations  of  mass  centers  and 
orientations  of  all  molecules  in  both  experimental  and 
predicted  unit  cells  are  calculated  and  compared. 
Orientations  of  the  molecules  are  described  in  terms  of  the 
Euler  angles  used  to  rotate  the  Cartesian  coordinate  axes 
onto  those  of  the  principal  axes  of  each  molecule. 
Finally,  a  visual  comparison  of  the  superimposed  unit  cells 
is  performed  to  confirm  that  the  MOLPAK  unit  cell 
resembles  the  experimental  cell. 

5.  RESULTS 

Several  types  of  data  are  extracted  from  the  ab  initio 
MP  results:  First,  we  determine  the  lattice  energies  of  the 
global  minimum  (GM)  and  that  of  the  structure  which 
represents  the  match  of  the  corresponding  experimental 
configuration.  This  last  type  of  structure  will  be  denoted 
in  the  following  as  “the  match”.  In  the  majority  of  cases 
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the  global  minimum  coincides  with  the  match  to 
experiment  and  as  a  result  the  two  energies  are  identical. 
We  also  calculated  the  energy  gap  between  the  global 
minimum  and  the  next  local  minimum  with  energy 
closest  to  the  GM  determined  from  the  entire  set  of 
hypothetical  crystal  structures  of  a  particular  system.  For 
each  particular  case  the  energetic  range  sampled  by  the 
corresponding  hypothetical  crystal  structures  was 
determined.  We  also  determined  the  percentage  of 
crystals  with  energies  smaller  than  that  of  the  match 
structure  and  with  energies  within  3  kcal/mol  from  the 
GM.  Finally,  we  determined  the  percentage  of  crystals 
with  energies  within  3  kcal/mol  of  GM  and  with 
densities  greater  than  the  one  determined  for  the 
structure  at  the  global  minimum.  All  of  this  information 
is  illustrated  in  Figure  1. 


Crystal  Index 

Figure  1  (a)  The  energy  range  of  the  hypothetical 

crystals  generated  for  each  system;  (b)  The  energy  gap 
AE  between  the  global  minimum  (GM)  and  the  local 
minimum  with  energy  closest  to  the  GM;  (c)  Percentages 
of  systems  that  are  within  3  kcal/mol  of  the  GM;  (d) 
Percentages  of  systems  that  have  lattice  energies  lower 
than  the  the  structure  that  matches  the  experiment;  (e) 
Percentages  of  systems  with  densities  higher  than  that  of 
the  GM.  The  crystal  indices  for  each  frame  correspond  to 
those  assigned  to  the  various  chemical  systems  listed  in 
Table  2S  of  Sorescu  and  Rice  (Sorescu  and  Rice,  2004). 
The  energy  units  are  kcal/mol. 


The  ranges  of  lattice  energies  relative  to  the  GM  for 
the  various  hypothetical  crystals  are  represented  in  Figure 
1(a).  It  can  be  seen  that  there  is  a  relatively  large  spread  of 
the  energetic  ranges  with  values  between  4.6  kcal/mol  to 
nearly  40  kcal/mol.  It  has  been  recommended  that 
structures  with  lattice  energies  that  are  more  than  3 
kcal/mol  higher  than  that  of  the  GM  should  not  be 
considered  as  probable  structures  (Motherwell  et  al., 
2002);  therefore,  for  those  systems  in  which  the  energy 
ranges  are  quite  large,  the  higher-energy  crystals  can  be 
discarded  from  consideration.  Furthermore,  it  has  been 
proposed  that  for  those  cases  in  which  the  lattice  energy  of 
the  match  to  experiment  is  the  GM,  a  significant  energy 
gap  between  the  GM  and  all  other  crystals  indicates  “a 
readily  predicted  robust  crystal  structure”  (Beyer  et  ah, 
2001).  The  corresponding  energy  gaps  obtained  in  our 
simulations  are  presented  in  Figure  1(b).  As  indicated  by 
this  figure  the  gap  energies  have  values  below  2  kcal/mol 
for  the  great  majority  of  crystals  investigated  indicating  the 
presence  of  closely  spaced  local  minima  in  the  region  of 
the  GM.  Also,  as  evident  Figure  1(b),  there  are  several 
cases  in  which  the  hypothetical  crystals  have  lattice 
energies  that  are  almost  equal  to  the  GM,  even  though 
upon  inspection,  the  crystal  parameters  were  often  quite 
dissimilar  (e.g.  different  space  groups,  different 
orientations  of  the  molecules  within  the  unit  cell). 

A  similar  finding  emerges  from  the  analysis  of  the 
data  illustrated  in  Figure  1(c)  where  the  percent  of 
hypothetical  crystals  with  lattice  energies  within  3 
kcal/mol  of  the  global  minimum  is  represented.  Indeed,  it 
can  be  observed  that  for  many  systems  a  large  number  of 
the  hypothetical  crystals  are  distributed  in  a  relative 
narrow  energy  window  above  the  GM. 

As  shown  in  Figure  1(d)  where  the  percentages  of 
crystals  that  are  lower  in  energy  than  the  matches  are 
given,  the  hypothetical  crystal  selected  as  the  match  to  the 
experimental  structure  is,  in  most  cases,  among  the  lowest- 
energy  crystals  generated  in  the  MOLPAK/WMIN 
procedure. 

We  have  investigated  whether  crystalline  polymorphs 
with  higher  densities  have  correspondingly  lower  lattice 
energies;  the  corresponding  results  are  presented  in  Figure 
1(e).  In  this  figure  we  represent  the  percentage  of 
hypothetical  crystals  with  densities  greater  than  that  of  the 
structure  corresponding  to  the  GM.  As  indicated  in  Figure 
1(e),  we  observe  that  in  majority  of  analyzed  cases, 
systems  having  the  highest  density  correspond  to  the  GM 
while  a  relatively  small  number  of  crystals  have  densities 
in  excess  of  that  of  GM. 

One  of  the  systems  having  a  GM  with  structural 
parameters  in  good  agreement  with  experiment  and  with  a 
significant  energy  gap  between  the  GM  and  all  the  other 
hypothetical  crystals  is  CTMTNA  (RDX, 
Cyclotrimethylene-trinitramine).  In  this  case  there  is  also 
a  moderate  range  of  spread  of  energies  (12.75  kcal/mol) 
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for  various  hypothetical  structures.  The  predicted 
densities  and  the  corresponding  lattice  energies  for  all 
local  energy  minima  identified  are  shown  in  Figure  2.  In 
this  case,  the  difference  in  energy  between  the  GM  and 
all  the  other  crystals  is  0.87  kcal/mol.  The  next  local 
minimum  above  the  GM  has  the  same  space  group 
symmetry  as  the  experimental  structure  and  has 
crystallographic  parameters  that  deviate  from  experiment 
by  only  a  few  percent.  However,  the  rms  rigid  body 
rotational  and  translational  displacements  for  this 
structure  are  large,  and  the  pictorial  representations  of 
the  hypothetical  crystal  superimposed  onto  the 
experimental  unit  cell  confirmed  that  this  higher-energy 
crystal  could  not  be  considered  as  a  good  match  for  this 
system. 


1.2  1.3  1.4  1.5  1.6  1.7  1.8 

Density  (g/cc) 

Figure  2  Lattice  energies  versus  densities  of 
hypothetical  crystals  of  METNAM22  and  RDX 
generated  by  ab  initio  crystal  prediction  methods. 
Crystalline  space  groups  in  which  minimizations  were 
performed  are  denoted  by  the  symbols  defined  in  the 
legend. 

For  the  RDX  system,  only  4.6%  of  the  crystals  have 
energies  within  3  kcal/mol  of  the  GM,  and  as  a  result  of 
lattice  minimization  48%  of  these  (namely  2.2%  from  the 
total  number)  have  converged  to  the  same  energy 
minimum  (the  GM)  after  minimization.  For  this  crystal 
it  has  also  been  found  that  the  structure  corresponding  to 
the  GM  has  also  the  highest  density  among  the  series  of 
hypothetical  crystals.  In  contradistinction  to  the  RDX 
case  we  have  also  found  several  situations  where  the 
lattice  energies  of  the  hypothetical  crystals  are  almost 
equal  to  the  GM  despite  significant  differences  in  the 
symmetries  and  crystallographic  parameters  of  these 
structures.  Such  a  case  is  the  METNAM22  crystal 
(Hexadeutero-N,N-dimethyl-nitramine,),  also  illustrated 
in  Figure  2.  In  this  case  the  range  of  lattice  energies  of 
the  predicted  hypothetical  crystals  is  4.5  kcal/mol  while 
99%  of  the  hypothetical  crystals  have  lattice  energies 
within  3  kcal/mol  of  the  global  minimum  (GM).  The 
GM  for  METNAM22  is  identified  as  the  Match  to  the 
experimental  crystal  having  the  monoclinic  P2!/c 


symmetry.  The  local  minimum  with  energy  closest  to  the 
GM  identified  by  MOLPAK/WMIN  procedure  has  space 
group  symmetry  PI  and  is  separated  by  only  0.02 
kcal/mol  from  the  GM.  Symmetry-constrained  normal 
mode  analyses  confirmed  that  these  energy  minima  are 
different. 


5.1.  Identification  of  Matches  to  Experimental 
Crystals 

In  most  cases,  the  best  match  of  MOLPAK  predictions 
to  experiment  is  easily  determined  and  the  superposition  of 
the  MOLPAK  match  to  experiment  is  so  good  that 
visually,  the  molecules  are  indistinguishable.  Such  a  case 
is,  for  example,  the  BECJEY  (2,3 -Dimethyl-2,3 - 
dinitrobutane)  crystal  represented  in  Figure  3a.  However, 
there  were  a  few  cases  in  which  the  molecular  orientations 
and  mass  positions  within  the  cell  had  larger  differences 
from  experiment.  Such  an  example  is  illustrated  in  Figure 
3b  for  the  case  of  the  DNEDAM  crystal. 


Figure  3  Pictorial  views  of  the  unit  cells  structures 
predicted  by  MOLPAK  superimposed  onto  the 
experimental  unit  cells  for  the  (a)  BECJEY  and  (b) 
DNEDAM  crystals.  The  cell  vectors  correspond  to  the 
experimental  values. 

The  MOLPAK/WMIN  calculations  generated  crystals 
that  matched  the  experimental  structure  for  148  of  the  174 
systems  studied  (85%).  Eight  of  the  148  matches  were 
found  only  after  performing  a  second  set  of  calculations  in 
which  the  test  molecule  is  reoriented  about  each  of  the 
Eulerian  axes  in  increments  of  5°  during  the  MOLPAK 
portion  of  the  calculations.  75%  (1 1 1)  of  the  148  matches 
have  the  lowest  lattice  energy  of  all  candidate  crystals 
generated  in  the  MOLPAK/WMIN  procedure. 

Another  point  we  have  investigated  in  this  study  is 
whether  the  aforementioned  statistics  on  the  capability  of 
the  MOLPAK/WMIN  procedure  to  identify  crystal 
structures  consistent  with  experiment  are  dependent  on  the 
chemical  classes  of  CHNO  crystals  included  in  this  study, 
i.e.  nitramines,  nitroaliphatics,  nitroaromatics  and  nitrate 
esters.  In  performing  this  analysis  we  need  to  point  out 
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that  several  systems  have  been  included  in  more  than  one 
category  since  many  of  them  have  multiple  NO2  groups 
that  are  attached  to  different  atom  types.  For  example, 
the  ENPROP  (Ethyl  3-nitrato-2-nitro-3-(4-nitrophenyl)- 
propionate)  crystal  has  N02  groups  attached  to  oxygen, 
an  aliphatic  carbon  and  an  aromatic  carbon.  Therefore, 
in  the  subsequent  analyses,  ENPROP  is  categorized  as  a 
nitrate  ester,  a  nitroaliphatic  and  a  nitroaromatic  system. 

Nitramines:  This  category  contains  77  crystals  that 
were  subjected  to  ab  initio  crystal  prediction.  Of  these, 
84%  (65)  were  identified  as  a  match  to  experiment  by  the 
MOLPAK/WMIN  procedure.  78%  (51)  of  the  crystals 
identified  as  matches  to  experiment  correspond  to  the 
global  energy  minimum. 

Nitroaliphatic:  85%  (56)  of  the  66  candidates  in 
this  category  were  identified  as  matches  to  the 
experimental  systems  by  the  MOLPAK/WMIN 
procedure;  of  these,  73%  (41)  correspond  to  the  global 
energy  minimum. 

Nitroaromatic:  27  out  of  29  of  the  systems  in  this 
category  (93%)  were  identified  as  matches  to 
experiment.  Of  these,  78%  (21)  correspond  to  the  global 
energy  minimum 

Nitrate  Esters:  Of  the  33  crystals  in  this  category, 
85%  (28)  were  identified  as  matches  to  the  experimental 
crystals;  79%  of  these  (22  out  of  28)  correspond  to  the 
global  energy  minimum. 

Separate  analyses  were  not  performed  for  systems 
that  could  be  categorized  as  furoxans  or  systems 
containing  azido  or  nitroso  groups,  due  to  the  small 
numbers  of  these  types  of  systems  studied.  The  above 
results  indicate  that  in  terms  of  identification  of  the 
experimental  structure,  the  method  and  SRT  model 
behave  fairly  consistently  across  all  categories 
investigated  with  a  success  rate  of  about  85%  for  the 
nitramine,  nitroaliphatic  and  nitrate  ester  categories  and 
an  even  slightly  higher  success  rate  for  the  nitroaromatic 
class  of  93%.  This  consistency  is  also  observed  with 
regard  to  the  percentages  of  the  identified  matched 
systems  corresponding  to  the  global  minimum.  In  this 
case  the  success  rate  ranges  from  73%  to  79%. 

For  the  37  matches  that  did  not  have  the  lowest 
lattice  energy  among  the  possible  candidates,  all  but  six 
had  lattice  energies  within  1  kcal/mol  of  the  global 
energy  minimum.  The  largest  difference  in  energy 
between  the  match  to  experiment  and  the  structure  with 
lowest  lattice  energy  is  for  PERYTN01  (Pentaerythritol 
tetranitrate,  2.86  kcal/mol).  This  particular  result  could 
be  indicative  of  deficiencies  in  the  SRT  force  field  for 
this  system.  However,  for  this  system,  no  R  factor  was 
recorded;  so  it  is  also  possible  that  the  structure  is  not 
well  characterized  at  the  experimental  conditions. 

An  additional  analysis  was  performed  to  determine 
whether  failures  of  the  search  method  clustered  into 


particular  space  groups.  For  the  most  part,  no  clustering  of 
failures  occurred  in  any  particular  space  group.  For  each 
individual  space  group  analyzed  by  MOLPAK  we 
determined  that  80  to  100%  of  the  systems  with  that  space 
group  produced  matches  to  the  experimental  structures. 
An  exception  from  this  statistics  was  seen  for  the  case  of 
Pca21  crystal  symmetry  which  had  a  success  rate  of  50%. 
However,  in  this  case  there  were  only  four  systems  with 
this  space  group  symmetry. 

5.2.  Structural  Predictions. 

The  MOLPAK/WMIN  predictions  of  the  densities  for 
crystals  identified  as  best  matches  to  experiment  are 
illustrated  in  Figures  4. 
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Figure  4  Calculated  percentage  errors  between  the 
predicted  and  experimental  values  of  the  densities  for  (a) 
the  entire  set  of  crystals;  (b)  nitramine  crystals;  (c) 
nitroaliphatic  crystals;  (d)  nitroaromatic  crystals;  and  (e) 
nitrate  ester  crystals.  The  crystal  indices  for  each  frame 
correspond  to  those  assigned  to  the  various  chemical 
systems  listed  in  Table  3S  of  Rice  and  Sorescu  (Rice  and 
Sorescu,  2004). 

Figure  4(a)  shows  the  percent  deviation  in  density  of 
the  MOLPAK  predictions  compared  to  experimental 
values  for  the  entire  set  of  molecules.  Overall,  the 
MOLPAK  densities  are  overestimated  on  average  by  2.8%, 
and  have  a  rms  deviation  of  3.9%.  Also,  in  the  great 
majority  of  cases,  namely  131  out  of  148,  the  predicted 
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crystal  densities  were  found  to  be  larger  than  the 
corresponding  experimental  values.  This  is  not 
surprising,  since  the  MOLPAK/WMIN  predictions 
reflect  simple  energy  minimizations  that  do  not  include 
thermal  effects,  and  the  majority  of  the  crystal  structural 
parameters  were  determined  at  room  temperature. 
Earlier  NPT-MD  simulations  for  a  few  of  the  systems 
studied  herein  predicted  expansion  of  the  lattice  (i.e.  a 
decrease  in  density)  when  thermal  effects  are  included 
(Sorescu  et  al.,  1997,  1998a,  1998b,  1999a).  Therefore, 
the  MOLPAK/WMIN  results  should  be  considered  as 
upper  limits  to  predictions  of  density,  assuming  that  no 
unusual  thermal  expansion  behavior  exists. 

The  percent  differences  in  density  between 
predictions  and  experiment  for  the  four  categories  of 
CHNO  crystals  analyzed  in  this  study  are  illustrated  in 
Fig.  4(b)-(e).  The  nitroaromatic  category  has  the 
smallest  rms  deviation  of  predictions  of  density  from 
experiment  (2.9%),  and  the  nitrate  ester  class  had  the 
largest  (4.7%).  The  rms  deviations  of  the  predictions  of 
density  from  experimental  values  for  the  nitramine  and 
nitroaliphatic  classes  are  3.6  and  4.0%,  respectively. 
Average  percent  deviations  of  density  between  the 
MOLPAK/WMIN  predictions  and  experiment  are  2.3, 

3.1,  2.3  and  4.1%  for  the  nitramine,  nitroaliphatic, 
nitroaromatic,  and  nitrate  ester  classes,  respectively. 

The  results  indicate  that  for  the  majority  of  the 
crystals,  the  predicted  lattice  dimensions  are  smaller  than 
experiment.  The  average  percent  differences  for  the  cell 
edge  lengths  a,  b  and  c  are  all  -1%,  while  the 
corresponding  rms  deviations  of  percent  differences  are 

2.2,  2.2  and  2.1%,  respectively.  Since  the  symmetry  of 
the  crystals  is  constrained  during  the  MOLPAK/WMIN 
calculations,  orthorhombic  and  monoclinic  systems  have 
some  or  all  cell  angles  fixed  at  90°  and  were  not  allowed 
to  vary  during  the  packing/minimization  procedure. 
Thus,  these  angles  were  not  included  in  the  distribution 
of  absolute  differences  in  predictions  of  cell  angles  from 
experimental  values.  The  total  number  of  non-fixed  cell 
angles  used  in  the  distribution  is  122.  Among  these 
about  88%  of  cell  angles  deviate  from  experiment  by  no 
more  than  2°. 

Regarding  molecular  orientation,  the  rms  rigid-body 
rotational  displacement  and  translational  displacement 
after  minimization  for  the  entire  set  of  crystals  are  1.9° 
and  0.07  A,  respectively.  The  rms  translational 
displacement  for  each  of  the  nitramine,  nitroaliphatic, 
nitroaromatic  and  nitrate  ester  classes  is  0.07  A.  The 
corresponding  rms  rigid-body  rotational  displacement  for 
same  four  categories,  are  2.1,  1.8,  2.0  and  1.6°, 
respectively. 

5.3,  Problem  cases 


We  have  indicated  in  previous  sections  that  the 
MOLPAK/WMIN  procedure  in  conjunction  with  the  SRT 
potential  was  able  to  match  the  experimental  structure  for 
148  out  of  174  systems  studied.  In  this  section  we  focus 
our  efforts  on  possible  sources  of  failure  for  the  remaining 
26  crystals.  Particularly,  in  an  attempt  to  determine 
whether  the  failure  to  identify  the  26  of  the  174  crystals 
was  due  to  inadequacies  in  the  SRT  intermolecular 
interaction  potential  or  in  the  search  method,  we  performed 
geometry  optimizations  of  the  26  crystals  using  the 
experimental  crystals  as  initial  structures  in  the  MP 
calculations.  Each  calculation  resulted  in  convergence  to 
a  local  energy  minimum  and  the  corresponding  crystal 
structure  was  superimposed  on  the  experimental  structure 
in  the  manner  described  in  Section  4.2.  Although  there  are 
notable  disagreements  in  the  predicted  crystallographic 
parameters  and  molecular  orientations  with  experiment, 
the  predicted  crystals  resemble  the  experimental  systems. 
In  all  cases,  the  predicted  crystal  densities  were  found  to 
be  smaller  than  the  corresponding  experimental  values  and 
on  average  underestimated  the  experimental  values  by 
2.9%.  This  is  in  contrast  to  the  results  for  systems 
identified  through  the  MOLPAK  procedure;  in  those  cases, 
89%  of  the  systems  had  densities  that  were  larger  than 
experiment.  The  rms  deviation  of  the  density  for  these 
problem  case  crystals  is  3.6%.  Average  percent 
differences  for  the  cell  edge  lengths  a,  b  and  c  are  -2.8,  - 
2.5,  and  -2.8%,  respectively  while  the  corresponding  rms 
deviations  of  percent  differences  are  4.3,  3.4  and  5.0%, 
respectively.  The  rms  rigid-body  rotational  and 
translational  displacements  after  minimization  for  this  set 
of  crystals  are  3.6°  and  0.13  A.  Overall,  these  results 
indicate  that  the  SRT  potential  did  a  poorer  job  of 
representing  these  systems.  Two  systems  were  particularly 
poorly  described  by  the  SRT  potential:  NOGUNA01  (N- 
Methyl-N-nitroso-N’-nitroguanidine)  and  JIHVEB  [(E)- 
Acetonitrolic  acid].  For  example,  we  have  seen  deviation 
of  the  cell  edge  lengths  as  high  as  19.7%  for  NOGUNA01 
and  10.76%  for  JIHVEB. 

Beside  deficiencies  of  the  potential  used  in  simulation 
of  these  crystals  we  note  that  another  reason  why  some  of 
the  searches  might  have  been  unsuccessful  is  due  to  the 
existence  of  several  local  minima  on  the  potential  energy 
surface  landscape.  As  a  result,  optimization  of  one  of  the 
candidate  crystals  generated  during  the  search  procedure 
might  be  prevented  from  reaching  the  GM  or  the  local 
minimum  corresponding  to  the  experimental  crystal  due  to 
convergence  to  one  of  the  other  local  minima.  In  such  an 
instance,  the  identification  of  a  local  minimum 
corresponding  to  the  experimental  structure  is  precluded. 
Beyer,  Lewis  and  Price  (Beyer  et  al.,  2001)  have  suggested 
that  the  search  procedure  is  flawed  for  those  cases  in 
which  a  local  minimum  with  structural  parameters  close  to 
experiment  exists  but  which  is  not  found  during  the  search 
procedure  (Beyer  et  al.,  2001;  Lommerse  et  al,  2000). 
Our  results  have  clearly  demonstrated  that  for  the  problem 
cases,  there  are  several  local  minima  with  structural 
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parameters  significantly  different  from  those  of  the 
experimental  structure.  Since  the  direct  minimization  of 
the  experimental  structure  for  each  of  the  problem  cases 
resulted  in  convergence  to  a  local  minimum  with 
structural  parameters  resembling  experiment,  the  failure 
to  identify  these  minima  in  the  MOLPAK/WMIN 
procedure  suggests  that  the  search  procedure  is  not  fully 
adequate  for  these  cases. 


6.  SUMMARY  AND  CONCLUSIONS 

A  total  of  174  CHNO  crystals  for  which 
experimental  crystallographic  information  is  available 
were  subjected  to  ab  initio  crystal  prediction  using  the 
MOLPAK/WMIN  suite  of  software  (Holden  et  al.,  1993; 
Busing,  1981)  and  the  Sorescu-Rice-Thompson  (SRT) 
interaction  potential  for  CHNO  crystals  (Sorescu  et  al., 
1997).  The  systems  under  study  consist  of  nitramine, 
nitroaliphatic,  nitroaromatic  and  nitrate  ester  molecules 
and  included  a  variety  of  acyclic,  monocyclic  and 
polycyclic/caged  species.  The  calculations  produced 
148  crystals  whose  crystallographic  parameters  and 
molecular  configurations  matched  those  of  the 
experimental  counterpart.  Of  these,  75%  corresponded 
to  the  global  energy  minimum  on  the  potential  energy 
surface.  Predicted  structural  parameters  of  the  matches 
were  in  good  agreement  with  experimental  values; 
predicted  densities  were  higher  than  experiment  by  less 
than  3%  on  average,  and  cell  edge  lengths  deviated  from 
experiment  on  average  by  no  more  than  -1%.  Rigid  body 
rotational  and  translational  displacements  had  rms 
deviations  from  experiment  of  1.9°  and  0.07  A, 
respectively. 

For  each  of  the  chemical  CHNO  classes  studied  in 
this  work,  namely  nitramines,  nitroaliphatic, 
nitroaromatic  and  nitrate  esters,  we  found  similar 
statistics  regarding  identification  of  crystal  and  structural 
parameters.  A  slightly  higher  percent  deviation  of 
density  from  experiment  was  found  for  the  nitrate  ester 
class.  The  uniformity  of  the  statistics  for  the  four 
categories  of  CHNO  crystals  that  were  identified  in  the 
MOLPAK/WMIN  procedure  indicates  that  the  SRT 
interaction  potential  is  transferable  among  these  classes 
of  compounds.  Consequently,  this  assessment  provides  a 
level  of  confidence  in  the  SRT  model  when  applied  to 
these  classes  of  CHNO  compounds.  The  overall  good 
performance  of  the  MOLPAK/WMIN  procedure  and 
SRT  potential  in  identifying  the  experimental  systems 
and  producing  relatively  good  agreement  with  structural 
parameters  for  85%  of  the  systems  studied  demonstrates 
that  both  the  model  and  method  are  reasonable  to  use  in 
ab  initio  crystal  predictions  of  similar  CHNO  systems 
when  the  molecular  configuration  is  known. 
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